

<!DOCTYPE html>
<!--[if IE 8]><html class="no-js lt-ie9" lang="en" > <![endif]-->
<!--[if gt IE 8]><!--> <html class="no-js" lang="en" > <!--<![endif]-->
<head>
  <meta charset="utf-8">
  
  <meta name="viewport" content="width=device-width, initial-scale=1.0">
  
  <title>quippy.elasticity &mdash; quippy 7ea04c5+ documentation</title>
  

  
  
    <link rel="shortcut icon" href="../../_static/favicon.ico"/>
  
  
  

  
  <script type="text/javascript" src="../../_static/js/modernizr.min.js"></script>
  
    
      <script type="text/javascript" id="documentation_options" data-url_root="../../" src="../../_static/documentation_options.js"></script>
        <script type="text/javascript" src="../../_static/jquery.js"></script>
        <script type="text/javascript" src="../../_static/underscore.js"></script>
        <script type="text/javascript" src="../../_static/doctools.js"></script>
        <script type="text/javascript" src="../../_static/language_data.js"></script>
        <script async="async" type="text/javascript" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/latest.js?config=TeX-AMS-MML_HTMLorMML"></script>
        <script type="text/x-mathjax-config">MathJax.Hub.Config({"tex2jax": {"processClass": "math|output_area", "processEscapes": true, "ignoreClass": "document", "inlineMath": [["$", "$"], ["\\(", "\\)"]]}})</script>
    
    <script type="text/javascript" src="../../_static/js/theme.js"></script>

    

  
  <link rel="stylesheet" href="../../_static/css/theme.css" type="text/css" />
  <link rel="stylesheet" href="../../_static/pygments.css" type="text/css" />
  <link rel="stylesheet" href="../../_static/theme_overrides.css" type="text/css" />
    <link rel="index" title="Index" href="../../genindex.html" />
    <link rel="search" title="Search" href="../../search.html" /> 
</head>

<body class="wy-body-for-nav">

   
  <div class="wy-grid-for-nav">
    
    <nav data-toggle="wy-nav-shift" class="wy-nav-side">
      <div class="wy-side-scroll">
        <div class="wy-side-nav-search" >
          

          
            <a href="../../index.html" class="icon icon-home"> quippy
          

          
            
            <img src="../../_static/hybrid.png" class="logo" alt="Logo"/>
          
          </a>

          
            
            
              <div class="version">
                7ea04c5+
              </div>
            
          

          
<div role="search">
  <form id="rtd-search-form" class="wy-form" action="../../search.html" method="get">
    <input type="text" name="q" placeholder="Search docs" />
    <input type="hidden" name="check_keywords" value="yes" />
    <input type="hidden" name="area" value="default" />
  </form>
</div>

          
        </div>

        <div class="wy-menu wy-menu-vertical" data-spy="affix" role="navigation" aria-label="main navigation">
          
    
            
            
              
            
            
              <ul>
<li class="toctree-l1"><a class="reference internal" href="../../intro.html">Introduction to <code class="docutils literal notranslate"><span class="pre">quippy</span></code></a></li>
<li class="toctree-l1"><a class="reference internal" href="../../install.html">Installation of QUIP and quippy</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../Tutorials/index.html">Tutorials</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../quippy.html">Quippy library reference</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../visualisation.html">Visualisation Tools</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../fortran_wrapper.html">Appendix: wrapping Fortran 90 Code</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../teach_sparse.html">Command line options of the teach_sparse main program</a></li>
</ul>

            
          
    <p></p>
    <ul>
        <li><a href="../../genindex.html">Index</a></li>
        <li><a href="../../py-modindex.html">Module Index</a></li>
        <li><a href="../../search.html">Search Page</a></li>
        <li><a href="../index.html">Module Source Code</a></li>
    </ul>
  
        </div>
      </div>
    </nav>

    <section data-toggle="wy-nav-shift" class="wy-nav-content-wrap">

      
      <nav class="wy-nav-top" aria-label="top navigation">
        
          <i data-toggle="wy-nav-top" class="fa fa-bars"></i>
          <a href="../../index.html">quippy</a>
        
      </nav>


      <div class="wy-nav-content">
        
        <div class="rst-content">
        
          















<div role="navigation" aria-label="breadcrumbs navigation">

  <ul class="wy-breadcrumbs">
    
      <li><a href="../../index.html">Docs</a> &raquo;</li>
        
          <li><a href="../index.html">Module code</a> &raquo;</li>
        
      <li>quippy.elasticity</li>
    
    
      <li class="wy-breadcrumbs-aside">
        
      </li>
    
  </ul>

  
  <hr/>
</div>
          <div role="main" class="document" itemscope="itemscope" itemtype="http://schema.org/Article">
           <div itemprop="articleBody">
            
  <h1>Source code for quippy.elasticity</h1><div class="highlight"><pre>
<span></span><span class="c1"># HQ XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX</span>
<span class="c1"># HQ X</span>
<span class="c1"># HQ X   quippy: Python interface to QUIP atomistic simulation library</span>
<span class="c1"># HQ X</span>
<span class="c1"># HQ X   Copyright James Kermode 2010</span>
<span class="c1"># HQ X</span>
<span class="c1"># HQ X   These portions of the source code are released under the GNU General</span>
<span class="c1"># HQ X   Public License, version 2, http://www.gnu.org/copyleft/gpl.html</span>
<span class="c1"># HQ X</span>
<span class="c1"># HQ X   If you would like to license the source code under different terms,</span>
<span class="c1"># HQ X   please contact James Kermode, james.kermode@gmail.com</span>
<span class="c1"># HQ X</span>
<span class="c1"># HQ X   When using this software, please cite the following reference:</span>
<span class="c1"># HQ X</span>
<span class="c1"># HQ X   http://www.jrkermode.co.uk/quippy</span>
<span class="c1"># HQ X</span>
<span class="c1"># HQ XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX</span>

<span class="kn">from</span> <span class="nn">quippy</span> <span class="k">import</span> <span class="n">available_modules</span>
<span class="k">if</span> <span class="s1">&#39;scipy&#39;</span> <span class="ow">in</span> <span class="n">available_modules</span><span class="p">:</span>
    <span class="kn">from</span> <span class="nn">scipy</span> <span class="k">import</span> <span class="n">stats</span>
<span class="kn">from</span> <span class="nn">farray</span> <span class="k">import</span> <span class="o">*</span>
<span class="kn">from</span> <span class="nn">quippy.atoms</span> <span class="k">import</span> <span class="n">Atoms</span>
<span class="kn">from</span> <span class="nn">quippy</span> <span class="k">import</span> <span class="n">_elasticity</span>
<span class="kn">from</span> <span class="nn">quippy._elasticity</span> <span class="k">import</span> <span class="o">*</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="n">__all__</span> <span class="o">=</span> <span class="n">_elasticity</span><span class="o">.</span><span class="n">__all__</span> <span class="o">+</span> <span class="p">[</span><span class="s1">&#39;strain_matrix&#39;</span><span class="p">,</span> <span class="s1">&#39;stress_matrix&#39;</span><span class="p">,</span> <span class="s1">&#39;strain_vector&#39;</span><span class="p">,</span>
                                 <span class="s1">&#39;stress_vector&#39;</span><span class="p">,</span> <span class="s1">&#39;fit_elastic_constants&#39;</span><span class="p">,</span>
                                 <span class="s1">&#39;elastic_constants&#39;</span><span class="p">,</span> <span class="s1">&#39;atomic_strain&#39;</span><span class="p">,</span>
                                 <span class="s1">&#39;elastic_fields_fortran&#39;</span><span class="p">,</span> <span class="s1">&#39;elastic_fields&#39;</span><span class="p">,</span>
                                 <span class="s1">&#39;AtomResolvedStressField&#39;</span><span class="p">,</span>
                                 <span class="s1">&#39;transform_elasticity&#39;</span><span class="p">,</span> <span class="s1">&#39;rayleigh_wave_speed&#39;</span><span class="p">]</span>

<div class="viewcode-block" id="strain_matrix"><a class="viewcode-back" href="../../elasticity.html#quippy.elasticity.strain_matrix">[docs]</a><span class="k">def</span> <span class="nf">strain_matrix</span><span class="p">(</span><span class="n">strain_vector</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Form a 3x3 strain matrix from a 6 component vector in Voigt notation</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">e1</span><span class="p">,</span> <span class="n">e2</span><span class="p">,</span> <span class="n">e3</span><span class="p">,</span> <span class="n">e4</span><span class="p">,</span> <span class="n">e5</span><span class="p">,</span> <span class="n">e6</span> <span class="o">=</span> <span class="n">strain_vector</span>
    <span class="k">return</span> <span class="n">farray</span><span class="p">([[</span><span class="mf">1.0</span><span class="o">+</span><span class="n">e1</span><span class="p">,</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">e6</span><span class="p">,</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">e5</span><span class="p">],</span>
                   <span class="p">[</span><span class="mf">0.5</span><span class="o">*</span><span class="n">e6</span><span class="p">,</span> <span class="mf">1.0</span><span class="o">+</span><span class="n">e2</span><span class="p">,</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">e4</span><span class="p">],</span>
                   <span class="p">[</span><span class="mf">0.5</span><span class="o">*</span><span class="n">e5</span><span class="p">,</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">e4</span><span class="p">,</span> <span class="mf">1.0</span><span class="o">+</span><span class="n">e3</span><span class="p">]])</span></div>

<div class="viewcode-block" id="stress_matrix"><a class="viewcode-back" href="../../elasticity.html#quippy.elasticity.stress_matrix">[docs]</a><span class="k">def</span> <span class="nf">stress_matrix</span><span class="p">(</span><span class="n">stress_vector</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Form a 3x3 stress matrix from a 6 component vector in Voigt notation</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">s1</span><span class="p">,</span> <span class="n">s2</span><span class="p">,</span> <span class="n">s3</span><span class="p">,</span> <span class="n">s4</span><span class="p">,</span> <span class="n">s5</span><span class="p">,</span> <span class="n">s6</span> <span class="o">=</span> <span class="n">stress_vector</span>
    <span class="k">return</span> <span class="n">farray</span><span class="p">([[</span><span class="n">s1</span><span class="p">,</span> <span class="n">s6</span><span class="p">,</span> <span class="n">s5</span><span class="p">],</span>
                   <span class="p">[</span><span class="n">s6</span><span class="p">,</span> <span class="n">s2</span><span class="p">,</span> <span class="n">s4</span><span class="p">],</span>
                   <span class="p">[</span><span class="n">s5</span><span class="p">,</span> <span class="n">s4</span><span class="p">,</span> <span class="n">s3</span><span class="p">]])</span></div>

<div class="viewcode-block" id="strain_vector"><a class="viewcode-back" href="../../elasticity.html#quippy.elasticity.strain_vector">[docs]</a><span class="k">def</span> <span class="nf">strain_vector</span><span class="p">(</span><span class="n">strain_matrix</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Form a 6 component strain vector in Voight notation from a 3x3 matrix</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="k">return</span> <span class="n">farray</span><span class="p">([</span><span class="n">strain_matrix</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="mf">1.0</span><span class="p">,</span>
                   <span class="n">strain_matrix</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span> <span class="o">-</span> <span class="mf">1.0</span><span class="p">,</span>
                   <span class="n">strain_matrix</span><span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">]</span> <span class="o">-</span> <span class="mf">1.0</span><span class="p">,</span>
                   <span class="mf">2.0</span><span class="o">*</span><span class="n">strain_matrix</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="mi">3</span><span class="p">],</span>
                   <span class="mf">2.0</span><span class="o">*</span><span class="n">strain_matrix</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">3</span><span class="p">],</span>
                   <span class="mf">2.0</span><span class="o">*</span><span class="n">strain_matrix</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">]])</span></div>

<div class="viewcode-block" id="stress_vector"><a class="viewcode-back" href="../../elasticity.html#quippy.elasticity.stress_vector">[docs]</a><span class="k">def</span> <span class="nf">stress_vector</span><span class="p">(</span><span class="n">stress_matrix</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Form a 6 component stress vector in Voight notation from a 3x3 matrix</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="k">return</span> <span class="n">farray</span><span class="p">([</span><span class="n">stress_matrix</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">],</span>
                   <span class="n">stress_matrix</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">],</span>
                   <span class="n">stress_matrix</span><span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">],</span>
                   <span class="n">stress_matrix</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="mi">3</span><span class="p">],</span>
                   <span class="n">stress_matrix</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">3</span><span class="p">],</span>
                   <span class="n">stress_matrix</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">]])</span></div>



<span class="n">Cij_symmetry</span> <span class="o">=</span> <span class="p">{</span>
   <span class="s1">&#39;cubic&#39;</span><span class="p">:</span>           <span class="n">farray</span><span class="p">([[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">7</span><span class="p">,</span> <span class="mi">7</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">7</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">7</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">7</span><span class="p">,</span> <span class="mi">7</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">4</span><span class="p">]]),</span>

   <span class="s1">&#39;trigonal_high&#39;</span><span class="p">:</span>   <span class="n">farray</span><span class="p">([[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">7</span><span class="p">,</span> <span class="mi">8</span><span class="p">,</span> <span class="mi">9</span><span class="p">,</span> <span class="mi">10</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">7</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">8</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span><span class="o">-</span><span class="mi">9</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">8</span><span class="p">,</span> <span class="mi">8</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">9</span><span class="p">,</span> <span class="o">-</span><span class="mi">9</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">10</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">6</span><span class="p">]]),</span>

   <span class="s1">&#39;trigonal_low&#39;</span><span class="p">:</span>    <span class="n">farray</span><span class="p">([[</span><span class="mi">1</span><span class="p">,</span>  <span class="mi">7</span><span class="p">,</span>  <span class="mi">8</span><span class="p">,</span>  <span class="mi">9</span><span class="p">,</span>  <span class="mi">10</span><span class="p">,</span>  <span class="mi">0</span> <span class="p">],</span>
                              <span class="p">[</span><span class="mi">7</span><span class="p">,</span>  <span class="mi">1</span><span class="p">,</span>  <span class="mi">8</span><span class="p">,</span> <span class="o">-</span><span class="mi">9</span><span class="p">,</span> <span class="o">-</span><span class="mi">10</span><span class="p">,</span>  <span class="mi">0</span> <span class="p">],</span>
                              <span class="p">[</span><span class="mi">8</span><span class="p">,</span>  <span class="mi">8</span><span class="p">,</span>  <span class="mi">3</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>   <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span> <span class="p">],</span>
                              <span class="p">[</span><span class="mi">9</span><span class="p">,</span> <span class="o">-</span><span class="mi">9</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">4</span><span class="p">,</span>   <span class="mi">0</span><span class="p">,</span> <span class="o">-</span><span class="mi">10</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">10</span><span class="p">,</span><span class="o">-</span><span class="mi">10</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>   <span class="mi">4</span><span class="p">,</span>  <span class="mi">9</span> <span class="p">],</span>
                              <span class="p">[</span><span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span> <span class="o">-</span><span class="mi">10</span> <span class="p">,</span> <span class="mi">9</span><span class="p">,</span>  <span class="mi">6</span> <span class="p">]]),</span>

   <span class="s1">&#39;tetragonal_high&#39;</span><span class="p">:</span> <span class="n">farray</span><span class="p">([[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">7</span><span class="p">,</span> <span class="mi">8</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">7</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">8</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">8</span><span class="p">,</span> <span class="mi">8</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">6</span><span class="p">]]),</span>

   <span class="s1">&#39;tetragonal_low&#39;</span><span class="p">:</span>  <span class="n">farray</span><span class="p">([[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">7</span><span class="p">,</span> <span class="mi">8</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">11</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">7</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">8</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="o">-</span><span class="mi">11</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">8</span><span class="p">,</span> <span class="mi">8</span><span class="p">,</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">11</span><span class="p">,</span> <span class="o">-</span><span class="mi">11</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">6</span><span class="p">]]),</span>

   <span class="s1">&#39;orthorhombic&#39;</span><span class="p">:</span>    <span class="n">farray</span><span class="p">([[</span> <span class="mi">1</span><span class="p">,</span>  <span class="mi">7</span><span class="p">,</span>  <span class="mi">8</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span> <span class="mi">7</span><span class="p">,</span>  <span class="mi">2</span><span class="p">,</span> <span class="mi">12</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span> <span class="mi">8</span><span class="p">,</span> <span class="mi">12</span><span class="p">,</span>  <span class="mi">3</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span> <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">4</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span> <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">5</span><span class="p">,</span>  <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span> <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">6</span><span class="p">]]),</span>

   <span class="s1">&#39;monoclinic&#39;</span><span class="p">:</span>      <span class="n">farray</span><span class="p">([[</span> <span class="mi">1</span><span class="p">,</span>  <span class="mi">7</span><span class="p">,</span>  <span class="mi">8</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">10</span><span class="p">,</span>  <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span> <span class="mi">7</span><span class="p">,</span>  <span class="mi">2</span><span class="p">,</span> <span class="mi">12</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span> <span class="mi">14</span><span class="p">,</span>  <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span> <span class="mi">8</span><span class="p">,</span> <span class="mi">12</span><span class="p">,</span>  <span class="mi">3</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span> <span class="mi">17</span><span class="p">,</span>  <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span> <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">4</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">20</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">10</span><span class="p">,</span> <span class="mi">14</span><span class="p">,</span> <span class="mi">17</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">5</span><span class="p">,</span>  <span class="mi">0</span><span class="p">],</span>
                              <span class="p">[</span> <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span> <span class="mi">20</span><span class="p">,</span>  <span class="mi">0</span><span class="p">,</span>  <span class="mi">6</span><span class="p">]]),</span>

   <span class="s1">&#39;triclinic&#39;</span><span class="p">:</span>       <span class="n">farray</span><span class="p">([[</span> <span class="mi">1</span><span class="p">,</span>  <span class="mi">7</span><span class="p">,</span>  <span class="mi">8</span><span class="p">,</span>  <span class="mi">9</span><span class="p">,</span>  <span class="mi">10</span><span class="p">,</span> <span class="mi">11</span><span class="p">],</span>
                              <span class="p">[</span> <span class="mi">7</span><span class="p">,</span>  <span class="mi">2</span><span class="p">,</span> <span class="mi">12</span><span class="p">,</span>  <span class="mi">13</span><span class="p">,</span> <span class="mi">14</span><span class="p">,</span> <span class="mi">15</span><span class="p">],</span>
                              <span class="p">[</span> <span class="mi">8</span><span class="p">,</span> <span class="mi">12</span><span class="p">,</span>  <span class="mi">3</span><span class="p">,</span>  <span class="mi">16</span><span class="p">,</span> <span class="mi">17</span><span class="p">,</span> <span class="mi">18</span><span class="p">],</span>
                              <span class="p">[</span> <span class="mi">9</span><span class="p">,</span> <span class="mi">13</span><span class="p">,</span> <span class="mi">16</span><span class="p">,</span>  <span class="mi">4</span><span class="p">,</span>  <span class="mi">19</span><span class="p">,</span> <span class="mi">20</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">10</span><span class="p">,</span> <span class="mi">14</span><span class="p">,</span> <span class="mi">17</span><span class="p">,</span> <span class="mi">19</span><span class="p">,</span>  <span class="mi">5</span><span class="p">,</span>  <span class="mi">21</span><span class="p">],</span>
                              <span class="p">[</span><span class="mi">11</span><span class="p">,</span> <span class="mi">15</span><span class="p">,</span> <span class="mi">18</span><span class="p">,</span> <span class="mi">20</span><span class="p">,</span>  <span class="mi">21</span><span class="p">,</span> <span class="mi">6</span> <span class="p">]]),</span>
   <span class="p">}</span>

<div class="viewcode-block" id="fit_elastic_constants"><a class="viewcode-back" href="../../elasticity.html#quippy.elasticity.fit_elastic_constants">[docs]</a><span class="k">def</span> <span class="nf">fit_elastic_constants</span><span class="p">(</span><span class="n">configs</span><span class="p">,</span> <span class="n">symmetry</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">N_steps</span><span class="o">=</span><span class="mi">5</span><span class="p">,</span> <span class="n">verbose</span><span class="o">=</span><span class="kc">True</span><span class="p">,</span> <span class="n">graphics</span><span class="o">=</span><span class="kc">True</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    quippy.elasticity.fit_elastic_constants() deprecated, please use</span>
<span class="sd">      matscipy.elasticity.fit_elastic_constants() instead</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="k">raise</span> <span class="ne">NotImplementedError</span><span class="p">(</span><span class="s2">&quot;quippy.elasticity.fit_elastic_contants() deprecated, please use matscipy.elasticity.fit_elastic_constants() instead&quot;</span><span class="p">)</span>    </div>

<div class="viewcode-block" id="elastic_constants"><a class="viewcode-back" href="../../elasticity.html#quippy.elasticity.elastic_constants">[docs]</a><span class="k">def</span> <span class="nf">elastic_constants</span><span class="p">(</span><span class="n">pot</span><span class="p">,</span> <span class="n">at</span><span class="p">,</span> <span class="n">sym</span><span class="o">=</span><span class="s1">&#39;cubic&#39;</span><span class="p">,</span> <span class="n">relax</span><span class="o">=</span><span class="kc">True</span><span class="p">,</span> <span class="n">verbose</span><span class="o">=</span><span class="kc">True</span><span class="p">,</span> <span class="n">graphics</span><span class="o">=</span><span class="kc">True</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    quippy.elasticity.elastic_contants() deprecated, please use</span>
<span class="sd">       matscipy.elasticity.fit_elastic_constants() instead</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="k">raise</span> <span class="ne">NotImplementedError</span><span class="p">(</span><span class="s2">&quot;quippy.elasticity.elastic_contants() deprecated, please use matscipy.elasticity.fit_elastic_constants() instead&quot;</span><span class="p">)</span></div>


<div class="viewcode-block" id="atomic_strain"><a class="viewcode-back" href="../../elasticity.html#quippy.elasticity.atomic_strain">[docs]</a><span class="k">def</span> <span class="nf">atomic_strain</span><span class="p">(</span><span class="n">at</span><span class="p">,</span> <span class="n">r0</span><span class="p">,</span> <span class="n">crystal_factor</span><span class="o">=</span><span class="mf">1.0</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;Atomic strain as defined by JA Zimmerman in `Continuum and Atomistic Modeling of</span>
<span class="sd">    Dislocation Nucleation at Crystal Surface Ledges`, PhD Thesis, Stanford University (1999).&quot;&quot;&quot;</span>

    <span class="n">strain</span> <span class="o">=</span> <span class="n">fzeros</span><span class="p">((</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">,</span><span class="n">at</span><span class="o">.</span><span class="n">n</span><span class="p">))</span>

    <span class="k">for</span> <span class="n">l</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="n">at</span><span class="o">.</span><span class="n">n</span><span class="p">):</span>
        <span class="k">for</span> <span class="n">n</span> <span class="ow">in</span> <span class="n">at</span><span class="o">.</span><span class="n">neighbours</span><span class="p">[</span><span class="n">l</span><span class="p">]:</span>
            <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
                <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
                    <span class="n">strain</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">l</span><span class="p">]</span> <span class="o">+=</span> <span class="n">n</span><span class="o">.</span><span class="n">diff</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">*</span><span class="n">n</span><span class="o">.</span><span class="n">diff</span><span class="p">[</span><span class="n">j</span><span class="p">]</span><span class="o">/</span><span class="n">r0</span><span class="o">**</span><span class="mf">2.0</span>

    <span class="k">return</span> <span class="n">strain</span><span class="o">/</span><span class="n">crystal_factor</span></div>


<span class="n">elastic_fields_fortran</span> <span class="o">=</span> <span class="n">elastic_fields</span> <span class="c1"># Fortran elastic fields routine</span>

<div class="viewcode-block" id="elastic_fields"><a class="viewcode-back" href="../../elasticity.html#quippy.elasticity.elastic_fields">[docs]</a><span class="k">def</span> <span class="nf">elastic_fields</span><span class="p">(</span><span class="n">at</span><span class="p">,</span> <span class="n">a</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span>  <span class="n">bond_length</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">c</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">c_vector</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">cij</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span>
                   <span class="n">save_reference</span><span class="o">=</span><span class="kc">False</span><span class="p">,</span> <span class="n">use_reference</span><span class="o">=</span><span class="kc">False</span><span class="p">,</span> <span class="n">mask</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">interpolate</span><span class="o">=</span><span class="kc">False</span><span class="p">,</span>
                   <span class="n">cutoff_factor</span><span class="o">=</span><span class="mf">1.2</span><span class="p">,</span> <span class="n">system</span><span class="o">=</span><span class="s1">&#39;tetrahedric&#39;</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute atomistic strain field and linear elastic stress response.</span>

<span class="sd">    Stress and strain are stored in compressed Voigt notation::</span>

<span class="sd">       at.strain[:,i] = [e_xx,e_yy,e_zz,e_yz,e_xz,e_xy]</span>
<span class="sd">       at.stress[:,i] = [sig_xx, sig_yy, sig_zz, sig_yz, sig_xz, sig_xy]</span>

<span class="sd">    so that sig = dot(C, strain) in the appropriate reference frame.</span>

<span class="sd">    Four-fold coordinate atoms within `at` are used to define</span>
<span class="sd">    tetrahedra. The deformation of each tetrahedra is determined</span>
<span class="sd">    relative to the ideal structure, using `a` as the cubic lattice</span>
<span class="sd">    constant (related to bond length by a factor :math:`sqrt{3}/4`).</span>
<span class="sd">    This deformation is then split into a strain and a rotation</span>
<span class="sd">    using a Polar decomposition.</span>

<span class="sd">    If `save_reference` or `use_reference` are True then `at` must have</span>
<span class="sd">    a `primitive_index` integer property which is different for each</span>
<span class="sd">    atom in the primitive unit cell. `save_reference` causes the</span>
<span class="sd">    local strain and rotation for one atom of each primitive type</span>
<span class="sd">    to be saved as entries `at.params`. Conversely, `use_reference` uses</span>
<span class="sd">    this information and to undo the local strain and rotation.</span>

<span class="sd">    The strain is then transformed into the crystal reference frame</span>
<span class="sd">    (i.e. x=100, y=010, z=001) to calculate the stress using the `cij`</span>
<span class="sd">    matrix of elastic constants. Finally the resulting stress is</span>
<span class="sd">    transformed back into the sample frame.</span>

<span class="sd">    The stress and strain tensor fields are interpolated to give</span>
<span class="sd">    values for atoms which are not four-fold coordinated (for example</span>
<span class="sd">    oxygen atoms in silica).</span>

<span class="sd">    Eigenvalues and eigenvectors of the stress are stored in the</span>
<span class="sd">    properties `stress_evals`,`stress_evec1`, `stress_evec2` and</span>
<span class="sd">    `stress_evec3`, ordered decreasingly by eigenvalue so that the</span>
<span class="sd">    principal eigenvalue and eigenvector are `stress_evals[1,:]` and</span>
<span class="sd">    `stress_evec1[:,i]` respectively.</span>

<span class="sd">    &quot;&quot;&quot;</span>

    <span class="k">def</span> <span class="nf">compute_stress_eig</span><span class="p">(</span><span class="n">i</span><span class="p">):</span>
        <span class="n">D</span><span class="p">,</span> <span class="n">SigEvecs</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">eig</span><span class="p">(</span><span class="n">stress_matrix</span><span class="p">(</span><span class="n">at</span><span class="o">.</span><span class="n">stress</span><span class="p">[:,</span><span class="n">i</span><span class="p">]))</span>

        <span class="c1"># Order by descending size of eigenvalues</span>
        <span class="n">sorted_evals</span><span class="p">,</span> <span class="n">order</span> <span class="o">=</span> <span class="nb">zip</span><span class="p">(</span><span class="o">*</span><span class="nb">sorted</span><span class="p">(</span><span class="nb">zip</span><span class="p">(</span><span class="n">D</span><span class="p">,</span> <span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">3</span><span class="p">]),</span><span class="n">reverse</span><span class="o">=</span><span class="kc">True</span><span class="p">))</span>

        <span class="n">at</span><span class="o">.</span><span class="n">stress_eval</span><span class="p">[:,</span><span class="n">i</span><span class="p">]</span>  <span class="o">=</span> <span class="n">sorted_evals</span>
        <span class="n">at</span><span class="o">.</span><span class="n">stress_evec1</span><span class="p">[:,</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">SigEvecs</span><span class="p">[:,</span><span class="n">order</span><span class="p">[</span><span class="mi">0</span><span class="p">]]</span>
        <span class="n">at</span><span class="o">.</span><span class="n">stress_evec2</span><span class="p">[:,</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">SigEvecs</span><span class="p">[:,</span><span class="n">order</span><span class="p">[</span><span class="mi">1</span><span class="p">]]</span>
        <span class="n">at</span><span class="o">.</span><span class="n">stress_evec3</span><span class="p">[:,</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">SigEvecs</span><span class="p">[:,</span><span class="n">order</span><span class="p">[</span><span class="mi">2</span><span class="p">]]</span>

    <span class="k">if</span> <span class="nb">sum</span><span class="p">([</span><span class="n">a</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">,</span> <span class="n">bond_length</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">])</span> <span class="o">!=</span> <span class="mi">1</span><span class="p">:</span>
        <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s1">&#39;One of lattice constant or bond length must be given&#39;</span><span class="p">)</span>

    <span class="k">if</span> <span class="n">a</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
        <span class="n">a</span> <span class="o">=</span> <span class="n">bond_length</span><span class="o">*</span><span class="mi">4</span><span class="o">/</span><span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mf">3.</span><span class="p">)</span>

    <span class="c1"># Use hysteretic connect as we want nearest neighbour connectivity only</span>
    <span class="n">at</span><span class="o">.</span><span class="n">calc_connect_hysteretic</span><span class="p">(</span><span class="n">cutoff_factor</span><span class="p">,</span> <span class="n">cutoff_factor</span><span class="p">)</span>

    <span class="k">if</span> <span class="p">(</span><span class="n">save_reference</span> <span class="ow">or</span> <span class="n">use_reference</span><span class="p">)</span> <span class="ow">and</span> <span class="ow">not</span> <span class="n">at</span><span class="o">.</span><span class="n">has_property</span><span class="p">(</span><span class="s1">&#39;primitive_index&#39;</span><span class="p">):</span>
        <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s1">&#39;Property &quot;primitive_index&quot; missing from Atoms object&#39;</span><span class="p">)</span>

    <span class="c1"># Add various properties to store results</span>
    <span class="n">at</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;strain&#39;</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">,</span> <span class="n">n_cols</span><span class="o">=</span><span class="mi">6</span><span class="p">)</span>
    <span class="k">if</span> <span class="n">cij</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span><span class="p">:</span>
        <span class="n">at</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;stress&#39;</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">,</span> <span class="n">n_cols</span><span class="o">=</span><span class="mi">6</span><span class="p">)</span>
        <span class="n">at</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;stress_eval&#39;</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">,</span> <span class="n">n_cols</span><span class="o">=</span><span class="mi">3</span><span class="p">)</span>
        <span class="n">at</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;stress_evec1&#39;</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">,</span> <span class="n">n_cols</span><span class="o">=</span><span class="mi">3</span><span class="p">)</span>
        <span class="n">at</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;stress_evec2&#39;</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">,</span> <span class="n">n_cols</span><span class="o">=</span><span class="mi">3</span><span class="p">)</span>
        <span class="n">at</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;stress_evec3&#39;</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">,</span> <span class="n">n_cols</span><span class="o">=</span><span class="mi">3</span><span class="p">)</span>
        <span class="n">at</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;strain_energy_density&#39;</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">)</span>

    <span class="n">rotXYZ</span> <span class="o">=</span> <span class="n">fzeros</span><span class="p">((</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">))</span>
    <span class="n">rotXYZ</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1.0</span>
    <span class="n">rotXYZ</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1.0</span>
    <span class="n">rotXYZ</span><span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1.0</span>
    <span class="n">E</span> <span class="o">=</span> <span class="n">fidentity</span><span class="p">(</span><span class="mi">3</span><span class="p">)</span>

    <span class="n">n_at</span><span class="o">=</span><span class="mi">0</span>

    <span class="c1"># Consider first atoms with four neighbours</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="n">at</span><span class="o">.</span><span class="n">n</span><span class="p">):</span>
        <span class="k">if</span> <span class="n">mask</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span> <span class="ow">and</span> <span class="ow">not</span> <span class="n">mask</span><span class="p">[</span><span class="n">i</span><span class="p">]:</span> <span class="k">continue</span>

        <span class="n">neighb</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">neighbours</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>

        <span class="k">if</span> <span class="p">(</span><span class="n">system</span> <span class="o">==</span> <span class="s1">&#39;anatase&#39;</span> <span class="ow">and</span> <span class="n">at</span><span class="o">.</span><span class="n">z</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">==</span> <span class="mi">22</span><span class="p">):</span>
           <span class="nb">print</span> <span class="n">i</span><span class="p">,</span> <span class="nb">len</span><span class="p">(</span><span class="n">neighb</span><span class="p">)</span>

        <span class="k">if</span> <span class="p">(</span><span class="n">system</span> <span class="o">==</span> <span class="s1">&#39;tetrahedric&#39;</span> <span class="ow">and</span> <span class="nb">len</span><span class="p">(</span><span class="n">neighb</span><span class="p">)</span> <span class="o">==</span> <span class="mi">4</span><span class="p">)</span> <span class="ow">or</span> <span class="p">(</span><span class="n">system</span> <span class="o">==</span> <span class="s1">&#39;anatase&#39;</span> <span class="ow">and</span> <span class="nb">len</span><span class="p">(</span><span class="n">neighb</span><span class="p">)</span> <span class="o">==</span> <span class="mi">6</span><span class="p">):</span>

            <span class="k">if</span> <span class="p">(</span><span class="n">system</span> <span class="o">==</span> <span class="s1">&#39;tetrahedric&#39;</span> <span class="ow">and</span> <span class="nb">len</span><span class="p">(</span><span class="n">neighb</span><span class="p">)</span> <span class="o">==</span> <span class="mi">4</span><span class="p">):</span>

                <span class="c1">#print &#39;\nProcessing atom&#39;, i</span>

                <span class="c1"># Consider neighbours in order of their index within the primitive cell</span>
                <span class="k">if</span> <span class="nb">hasattr</span><span class="p">(</span><span class="n">at</span><span class="p">,</span> <span class="s1">&#39;primitive_index&#39;</span><span class="p">):</span>
                    <span class="p">(</span><span class="n">j1</span><span class="p">,</span><span class="n">i1</span><span class="p">),</span> <span class="p">(</span><span class="n">j2</span><span class="p">,</span><span class="n">i2</span><span class="p">),</span> <span class="p">(</span><span class="n">j3</span><span class="p">,</span><span class="n">i3</span><span class="p">),</span> <span class="p">(</span><span class="n">j4</span><span class="p">,</span><span class="n">i4</span><span class="p">)</span> <span class="o">=</span> <span class="nb">sorted</span><span class="p">((</span><span class="n">at</span><span class="o">.</span><span class="n">primitive_index</span><span class="p">[</span><span class="n">n</span><span class="o">.</span><span class="n">j</span><span class="p">],</span><span class="n">i</span><span class="p">)</span> <span class="k">for</span> <span class="n">i</span><span class="p">,</span><span class="n">n</span> <span class="ow">in</span> <span class="n">fenumerate</span><span class="p">(</span><span class="n">neighb</span><span class="p">))</span>
                <span class="k">else</span><span class="p">:</span>
                    <span class="p">(</span><span class="n">j1</span><span class="p">,</span><span class="n">i1</span><span class="p">),</span> <span class="p">(</span><span class="n">j2</span><span class="p">,</span><span class="n">i2</span><span class="p">),</span> <span class="p">(</span><span class="n">j3</span><span class="p">,</span><span class="n">i3</span><span class="p">),</span> <span class="p">(</span><span class="n">j4</span><span class="p">,</span><span class="n">i4</span><span class="p">)</span> <span class="o">=</span> <span class="nb">list</span><span class="p">((</span><span class="n">n</span><span class="o">.</span><span class="n">j</span><span class="p">,</span><span class="n">i</span><span class="p">)</span> <span class="k">for</span> <span class="n">i</span><span class="p">,</span><span class="n">n</span> <span class="ow">in</span> <span class="n">fenumerate</span><span class="p">(</span><span class="n">neighb</span><span class="p">))</span>

                <span class="c1"># Find cubic axes from neighbours</span>
                <span class="n">n1</span> <span class="o">=</span> <span class="n">neighb</span><span class="p">[</span><span class="n">i2</span><span class="p">]</span><span class="o">.</span><span class="n">diff</span> <span class="o">-</span> <span class="n">neighb</span><span class="p">[</span><span class="n">i1</span><span class="p">]</span><span class="o">.</span><span class="n">diff</span>
                <span class="n">n2</span> <span class="o">=</span> <span class="n">neighb</span><span class="p">[</span><span class="n">i3</span><span class="p">]</span><span class="o">.</span><span class="n">diff</span> <span class="o">-</span> <span class="n">neighb</span><span class="p">[</span><span class="n">i1</span><span class="p">]</span><span class="o">.</span><span class="n">diff</span>
                <span class="n">n3</span> <span class="o">=</span> <span class="n">neighb</span><span class="p">[</span><span class="n">i4</span><span class="p">]</span><span class="o">.</span><span class="n">diff</span> <span class="o">-</span> <span class="n">neighb</span><span class="p">[</span><span class="n">i1</span><span class="p">]</span><span class="o">.</span><span class="n">diff</span>
                <span class="c1">#print &#39;n1&#39;, n1.norm(), n1</span>
                <span class="c1">#print &#39;n2&#39;, n2.norm(), n2</span>
                <span class="c1">#print &#39;n3&#39;, n3.norm(), n3</span>

                <span class="n">E</span><span class="p">[:,</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">n1</span> <span class="o">+</span> <span class="n">n2</span> <span class="o">-</span> <span class="n">n3</span><span class="p">)</span><span class="o">/</span><span class="n">a</span>
                <span class="n">E</span><span class="p">[:,</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">n2</span> <span class="o">+</span> <span class="n">n3</span> <span class="o">-</span> <span class="n">n1</span><span class="p">)</span><span class="o">/</span><span class="n">a</span>
                <span class="n">E</span><span class="p">[:,</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">n3</span> <span class="o">+</span> <span class="n">n1</span> <span class="o">-</span> <span class="n">n2</span><span class="p">)</span><span class="o">/</span><span class="n">a</span>

                <span class="k">if</span> <span class="nb">all</span><span class="p">(</span><span class="n">E</span><span class="p">[:,</span><span class="mi">1</span><span class="p">]</span> <span class="o">==</span> <span class="mi">0</span><span class="p">)</span> <span class="ow">or</span> <span class="nb">all</span><span class="p">(</span><span class="n">E</span><span class="p">[:,</span><span class="mi">2</span><span class="p">]</span> <span class="o">==</span> <span class="mi">0</span><span class="p">)</span> <span class="ow">or</span> <span class="nb">all</span><span class="p">(</span><span class="n">E</span><span class="p">[:,</span><span class="mi">3</span><span class="p">])</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
                    <span class="k">continue</span>

            <span class="k">elif</span> <span class="p">(</span><span class="n">system</span> <span class="o">==</span> <span class="s1">&#39;anatase&#39;</span> <span class="ow">and</span> <span class="nb">len</span><span class="p">(</span><span class="n">neighb</span><span class="p">)</span> <span class="o">==</span> <span class="mi">6</span><span class="p">):</span>
                <span class="c1">#c_vector indicates the direction of the longest O-O bond (among the ones with Ti as a center of symmetry).</span>
                <span class="c1">#It is necessary to identify the crystal orientation. </span>

                <span class="n">n_at</span> <span class="o">=</span> <span class="n">n_at</span> <span class="o">+</span> <span class="mi">1</span>

                <span class="k">if</span><span class="p">(</span><span class="n">c_vector</span><span class="o">==</span><span class="kc">None</span><span class="p">):</span>
                    <span class="n">c_vector</span> <span class="o">=</span> <span class="n">fzeros</span><span class="p">(</span><span class="mi">3</span><span class="p">)</span>
                    <span class="n">c_vector</span><span class="p">[:]</span> <span class="o">=</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span>

                <span class="n">c_local</span> <span class="o">=</span> <span class="n">fzeros</span><span class="p">(</span><span class="mi">3</span><span class="p">)</span>
                <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
                    <span class="n">c_local</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">c_vector</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">*</span> <span class="p">(</span><span class="n">c</span><span class="o">/</span><span class="mi">2</span><span class="p">)</span>

                <span class="c1"># Consider neighbours in order of their index within the primitive cell</span>
                <span class="k">if</span> <span class="nb">hasattr</span><span class="p">(</span><span class="n">at</span><span class="p">,</span> <span class="s1">&#39;primitive_index&#39;</span><span class="p">):</span>
                    <span class="p">(</span><span class="n">j1</span><span class="p">,</span><span class="n">i1</span><span class="p">),</span> <span class="p">(</span><span class="n">j2</span><span class="p">,</span><span class="n">i2</span><span class="p">),</span> <span class="p">(</span><span class="n">j3</span><span class="p">,</span><span class="n">i3</span><span class="p">),</span> <span class="p">(</span><span class="n">j4</span><span class="p">,</span><span class="n">i4</span><span class="p">),</span> <span class="p">(</span><span class="n">j5</span><span class="p">,</span><span class="n">i5</span><span class="p">),</span> <span class="p">(</span><span class="n">j6</span><span class="p">,</span><span class="n">i6</span><span class="p">)</span> <span class="o">=</span> <span class="nb">sorted</span><span class="p">((</span><span class="n">at</span><span class="o">.</span><span class="n">primitive_index</span><span class="p">[</span><span class="n">n</span><span class="o">.</span><span class="n">j</span><span class="p">],</span><span class="n">i</span><span class="p">)</span> <span class="k">for</span> <span class="n">i</span><span class="p">,</span><span class="n">n</span> <span class="ow">in</span> <span class="n">fenumerate</span><span class="p">(</span><span class="n">neighb</span><span class="p">))</span>
                <span class="k">else</span><span class="p">:</span>
                    <span class="p">(</span><span class="n">j1</span><span class="p">,</span><span class="n">i1</span><span class="p">),</span> <span class="p">(</span><span class="n">j2</span><span class="p">,</span><span class="n">i2</span><span class="p">),</span> <span class="p">(</span><span class="n">j3</span><span class="p">,</span><span class="n">i3</span><span class="p">),</span> <span class="p">(</span><span class="n">j4</span><span class="p">,</span><span class="n">i4</span><span class="p">),</span> <span class="p">(</span><span class="n">j5</span><span class="p">,</span><span class="n">i5</span><span class="p">),</span> <span class="p">(</span><span class="n">j6</span><span class="p">,</span><span class="n">i6</span><span class="p">)</span> <span class="o">=</span> <span class="nb">list</span><span class="p">((</span><span class="n">n</span><span class="o">.</span><span class="n">j</span><span class="p">,</span><span class="n">i</span><span class="p">)</span> <span class="k">for</span> <span class="n">i</span><span class="p">,</span><span class="n">n</span> <span class="ow">in</span> <span class="n">fenumerate</span><span class="p">(</span><span class="n">neighb</span><span class="p">))</span>

                <span class="n">dd</span> <span class="o">=</span> <span class="n">fzeros</span><span class="p">(</span><span class="mi">6</span><span class="p">)</span>
                <span class="n">dd</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">norm</span><span class="p">(</span><span class="n">c_local</span> <span class="o">-</span> <span class="n">neighb</span><span class="p">[</span><span class="n">i1</span><span class="p">]</span><span class="o">.</span><span class="n">diff</span><span class="p">)</span>
                <span class="n">dd</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">norm</span><span class="p">(</span><span class="n">c_local</span> <span class="o">-</span> <span class="n">neighb</span><span class="p">[</span><span class="n">i2</span><span class="p">]</span><span class="o">.</span><span class="n">diff</span><span class="p">)</span>
                <span class="n">dd</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">norm</span><span class="p">(</span><span class="n">c_local</span> <span class="o">-</span> <span class="n">neighb</span><span class="p">[</span><span class="n">i3</span><span class="p">]</span><span class="o">.</span><span class="n">diff</span><span class="p">)</span>
                <span class="n">dd</span><span class="p">[</span><span class="mi">4</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">norm</span><span class="p">(</span><span class="n">c_local</span> <span class="o">-</span> <span class="n">neighb</span><span class="p">[</span><span class="n">i4</span><span class="p">]</span><span class="o">.</span><span class="n">diff</span><span class="p">)</span>
                <span class="n">dd</span><span class="p">[</span><span class="mi">5</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">norm</span><span class="p">(</span><span class="n">c_local</span> <span class="o">-</span> <span class="n">neighb</span><span class="p">[</span><span class="n">i5</span><span class="p">]</span><span class="o">.</span><span class="n">diff</span><span class="p">)</span>
                <span class="n">dd</span><span class="p">[</span><span class="mi">6</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">norm</span><span class="p">(</span><span class="n">c_local</span> <span class="o">-</span> <span class="n">neighb</span><span class="p">[</span><span class="n">i6</span><span class="p">]</span><span class="o">.</span><span class="n">diff</span><span class="p">)</span>
                <span class="c1">#print dd[5], c_local - neighb[i5].diff</span>
                <span class="n">ind</span> <span class="o">=</span>  <span class="n">np</span><span class="o">.</span><span class="n">argsort</span><span class="p">(</span><span class="n">dd</span><span class="p">)</span>
               
                <span class="n">dd2</span> <span class="o">=</span> <span class="n">fzeros</span><span class="p">(</span><span class="mi">3</span><span class="p">)</span>
                <span class="n">dd2</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">distance_min_image</span><span class="p">(</span> <span class="n">neighb</span><span class="p">[</span><span class="n">ind</span><span class="p">[</span><span class="mi">2</span><span class="p">]]</span><span class="o">.</span><span class="n">j</span> <span class="p">,</span> <span class="n">neighb</span><span class="p">[</span><span class="n">ind</span><span class="p">[</span><span class="mi">3</span><span class="p">]]</span><span class="o">.</span><span class="n">j</span> <span class="p">)</span>
                <span class="n">dd2</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">distance_min_image</span><span class="p">(</span> <span class="n">neighb</span><span class="p">[</span><span class="n">ind</span><span class="p">[</span><span class="mi">2</span><span class="p">]]</span><span class="o">.</span><span class="n">j</span> <span class="p">,</span> <span class="n">neighb</span><span class="p">[</span><span class="n">ind</span><span class="p">[</span><span class="mi">4</span><span class="p">]]</span><span class="o">.</span><span class="n">j</span> <span class="p">)</span>
                <span class="n">dd2</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">distance_min_image</span><span class="p">(</span> <span class="n">neighb</span><span class="p">[</span><span class="n">ind</span><span class="p">[</span><span class="mi">2</span><span class="p">]]</span><span class="o">.</span><span class="n">j</span> <span class="p">,</span> <span class="n">neighb</span><span class="p">[</span><span class="n">ind</span><span class="p">[</span><span class="mi">5</span><span class="p">]]</span><span class="o">.</span><span class="n">j</span> <span class="p">)</span>
                <span class="n">ind2</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">argsort</span><span class="p">(</span><span class="n">dd2</span><span class="p">)</span> 
                <span class="k">if</span>  <span class="n">ind2</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span> <span class="o">==</span> <span class="mi">2</span><span class="p">:</span>
                     <span class="n">a1</span> <span class="o">=</span> <span class="n">ind</span><span class="p">[</span><span class="mi">4</span><span class="p">]</span> 
                     <span class="n">a2</span> <span class="o">=</span> <span class="n">ind</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span> 
                     <span class="n">ind</span><span class="p">[</span><span class="mi">4</span><span class="p">]</span> <span class="o">=</span> <span class="n">a2</span> 
                     <span class="n">ind</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="n">a1</span>
                <span class="k">elif</span> <span class="n">ind2</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span> <span class="o">==</span> <span class="mi">3</span><span class="p">:</span>
                     <span class="n">a1</span> <span class="o">=</span> <span class="n">ind</span><span class="p">[</span><span class="mi">5</span><span class="p">]</span> 
                     <span class="n">a2</span> <span class="o">=</span> <span class="n">ind</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span> 
                     <span class="n">ind</span><span class="p">[</span><span class="mi">5</span><span class="p">]</span> <span class="o">=</span> <span class="n">a2</span> 
                     <span class="n">ind</span><span class="p">[</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="n">a1</span>

                <span class="n">n1</span> <span class="o">=</span> <span class="n">neighb</span><span class="p">[</span><span class="n">ind</span><span class="p">[</span><span class="mi">1</span><span class="p">]]</span><span class="o">.</span><span class="n">diff</span> <span class="o">-</span> <span class="n">neighb</span><span class="p">[</span><span class="n">ind</span><span class="p">[</span><span class="mi">6</span><span class="p">]]</span><span class="o">.</span><span class="n">diff</span>
                <span class="n">n2</span> <span class="o">=</span> <span class="n">neighb</span><span class="p">[</span><span class="n">ind</span><span class="p">[</span><span class="mi">2</span><span class="p">]]</span><span class="o">.</span><span class="n">diff</span> <span class="o">-</span> <span class="n">neighb</span><span class="p">[</span><span class="n">ind</span><span class="p">[</span><span class="mi">3</span><span class="p">]]</span><span class="o">.</span><span class="n">diff</span>
                <span class="n">n3</span> <span class="o">=</span> <span class="n">neighb</span><span class="p">[</span><span class="n">ind</span><span class="p">[</span><span class="mi">4</span><span class="p">]]</span><span class="o">.</span><span class="n">diff</span> <span class="o">-</span> <span class="n">neighb</span><span class="p">[</span><span class="n">ind</span><span class="p">[</span><span class="mi">5</span><span class="p">]]</span><span class="o">.</span><span class="n">diff</span>

                <span class="n">E</span><span class="p">[:,</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">n3</span><span class="o">/</span><span class="n">a</span>
                <span class="n">E</span><span class="p">[:,</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="n">n2</span><span class="o">/</span><span class="n">a</span>
                <span class="n">E</span><span class="p">[:,</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="n">n1</span><span class="o">/</span><span class="n">c</span>

            <span class="c1">#print &#39;E&#39;, E</span>

            <span class="c1"># Kill near zero elements</span>
            <span class="n">E</span><span class="p">[</span><span class="nb">abs</span><span class="p">(</span><span class="n">E</span><span class="p">)</span> <span class="o">&lt;</span> <span class="mf">1e-6</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.0</span>

            <span class="k">if</span> <span class="p">(</span><span class="n">E</span> <span class="o">&lt;</span> <span class="mf">0.0</span><span class="p">)</span><span class="o">.</span><span class="n">all</span><span class="p">():</span> <span class="n">E</span> <span class="o">=</span> <span class="o">-</span><span class="n">E</span>

            <span class="c1"># Find polar decomposition: E = S*R where S is symmetric,</span>
            <span class="c1"># and R is a rotation</span>
            <span class="c1">#</span>
            <span class="c1">#  EEt = E*E&#39;, EEt = VDV&#39; D diagonal, S = V D^1/2 V&#39;, R = S^-1*E</span>
            <span class="n">EEt</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">E</span><span class="p">,</span> <span class="n">E</span><span class="o">.</span><span class="n">T</span><span class="p">)</span>

            <span class="n">D</span><span class="p">,</span> <span class="n">V</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">eig</span><span class="p">(</span><span class="n">EEt</span><span class="p">)</span>

            <span class="n">S</span> <span class="o">=</span> <span class="n">farray</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">V</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">diag</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">D</span><span class="p">))),</span> <span class="n">V</span><span class="o">.</span><span class="n">T</span><span class="p">))</span>
            <span class="n">R</span> <span class="o">=</span> <span class="n">farray</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">V</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">diag</span><span class="p">(</span><span class="n">D</span><span class="o">**-</span><span class="mf">0.5</span><span class="p">)),</span> <span class="n">V</span><span class="o">.</span><span class="n">T</span><span class="p">),</span> <span class="n">E</span><span class="p">))</span>

            <span class="c1">#print &#39;S:&#39;, S</span>
            <span class="c1">#print &#39;R:&#39;, R</span>

            <span class="k">if</span> <span class="n">save_reference</span><span class="p">:</span>
                <span class="n">key</span> <span class="o">=</span> <span class="s1">&#39;strain_inv_</span><span class="si">%d</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">at</span><span class="o">.</span><span class="n">primitive_index</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
                <span class="k">if</span> <span class="n">key</span> <span class="ow">not</span> <span class="ow">in</span> <span class="n">at</span><span class="o">.</span><span class="n">params</span><span class="p">:</span>
                    <span class="n">at</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="n">key</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">inv</span><span class="p">(</span><span class="n">S</span><span class="p">)</span>

            <span class="k">if</span> <span class="n">use_reference</span><span class="p">:</span>
                <span class="n">S</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">S</span><span class="p">,</span> <span class="n">at</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;strain_inv_</span><span class="si">%d</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">at</span><span class="o">.</span><span class="n">primitive_index</span><span class="p">[</span><span class="n">i</span><span class="p">]])</span>

            <span class="c1">#print &#39;S after apply S0_inv:&#39;, S</span>

            <span class="c1"># Strain in rotated coordinate system</span>
            <span class="n">at</span><span class="o">.</span><span class="n">strain</span><span class="p">[:,</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">strain_vector</span><span class="p">(</span><span class="n">S</span><span class="p">)</span>

            <span class="c1"># Test for permutations - check which way x points</span>
            <span class="n">RtE</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">R</span><span class="o">.</span><span class="n">T</span><span class="p">,</span> <span class="n">E</span><span class="p">)</span>
            <span class="k">if</span> <span class="n">RtE</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">&gt;</span> <span class="n">RtE</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="ow">and</span> <span class="n">RtE</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">&gt;</span> <span class="n">RtE</span><span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">1</span><span class="p">]:</span>
                <span class="n">R</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">rotXYZ</span><span class="p">,</span> <span class="n">R</span><span class="p">)</span>
            <span class="k">elif</span> <span class="n">RtE</span><span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">&gt;</span> <span class="n">RtE</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="ow">and</span> <span class="n">RtE</span><span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">&gt;</span> <span class="n">RtE</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="mi">1</span><span class="p">]:</span>
                <span class="n">R</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">rotXYZ</span><span class="o">.</span><span class="n">T</span><span class="p">,</span> <span class="n">R</span><span class="p">)</span>

            <span class="k">if</span> <span class="n">save_reference</span><span class="p">:</span>
                <span class="n">key</span> <span class="o">=</span> <span class="s1">&#39;rotation_inv_</span><span class="si">%d</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">at</span><span class="o">.</span><span class="n">primitive_index</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
                <span class="k">if</span> <span class="n">key</span> <span class="ow">not</span> <span class="ow">in</span> <span class="n">at</span><span class="o">.</span><span class="n">params</span><span class="p">:</span>
                    <span class="n">at</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="n">key</span><span class="p">]</span> <span class="o">=</span> <span class="n">R</span><span class="o">.</span><span class="n">T</span>

            <span class="k">if</span> <span class="n">use_reference</span><span class="p">:</span>
                <span class="n">R</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">R</span><span class="p">,</span> <span class="n">at</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;rotation_inv_</span><span class="si">%d</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">at</span><span class="o">.</span><span class="n">primitive_index</span><span class="p">[</span><span class="n">i</span><span class="p">]])</span>

            <span class="c1">#print &#39;R after apply R0t:&#39;, R</span>

            <span class="c1"># Rotate to crystal coordinate system to apply Cij matrix</span>
            <span class="n">RtSR</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">R</span><span class="o">.</span><span class="n">T</span><span class="p">,</span> <span class="n">S</span><span class="p">),</span> <span class="n">R</span><span class="p">)</span>
            <span class="c1">#print &#39;RtSR:&#39;, RtSR</span>

            <span class="k">if</span> <span class="n">cij</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span><span class="p">:</span>
                <span class="n">sig</span> <span class="o">=</span> <span class="n">stress_matrix</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">cij</span><span class="p">,</span> <span class="n">strain_vector</span><span class="p">(</span><span class="n">RtSR</span><span class="p">)))</span>

                <span class="c1"># Rotate back to local coordinate system</span>
                <span class="n">RsigRt</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">R</span><span class="p">,</span> <span class="n">sig</span><span class="p">),</span> <span class="n">R</span><span class="o">.</span><span class="n">T</span><span class="p">)</span>

                <span class="c1"># Symmetrise stress tensor</span>
                <span class="n">RsigRt</span> <span class="o">=</span> <span class="p">(</span><span class="n">RsigRt</span> <span class="o">+</span> <span class="n">RsigRt</span><span class="o">.</span><span class="n">T</span><span class="p">)</span><span class="o">/</span><span class="mf">2.0</span>
                <span class="n">at</span><span class="o">.</span><span class="n">stress</span><span class="p">[:,</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">stress_vector</span><span class="p">(</span><span class="n">RsigRt</span><span class="p">)</span>

                <span class="n">at</span><span class="o">.</span><span class="n">strain_energy_density</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">at</span><span class="o">.</span><span class="n">strain</span><span class="p">[:,</span><span class="n">i</span><span class="p">],</span> <span class="n">at</span><span class="o">.</span><span class="n">stress</span><span class="p">[:,</span><span class="n">i</span><span class="p">])</span>
                <span class="c1">#print at.strain[:, i]</span>

                <span class="n">compute_stress_eig</span><span class="p">(</span><span class="n">i</span><span class="p">)</span>

    <span class="c1"># For atoms without 4 neighbours, interpolate stress and strain fields</span>
    <span class="k">if</span> <span class="n">interpolate</span><span class="p">:</span>
        <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="n">at</span><span class="o">.</span><span class="n">n</span><span class="p">):</span>
            <span class="k">if</span> <span class="n">mask</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span> <span class="ow">and</span> <span class="ow">not</span> <span class="n">mask</span><span class="p">[</span><span class="n">i</span><span class="p">]:</span> <span class="k">continue</span>

            <span class="n">neighb</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">neighbours</span><span class="p">[</span><span class="n">i</span><span class="p">]</span>
            <span class="k">if</span> <span class="n">at</span><span class="o">.</span><span class="n">z</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">!=</span> <span class="mi">8</span> <span class="ow">or</span>  <span class="nb">len</span><span class="p">(</span><span class="n">neighb</span><span class="p">)</span> <span class="o">==</span> <span class="mi">0</span> <span class="p">:</span> <span class="k">continue</span>

            <span class="n">at</span><span class="o">.</span><span class="n">strain</span><span class="p">[:,</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">strain</span><span class="p">[:,[</span><span class="n">n</span><span class="o">.</span><span class="n">j</span> <span class="k">for</span> <span class="n">n</span> <span class="ow">in</span> <span class="n">neighb</span><span class="p">]]</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">2</span><span class="p">)</span>
            <span class="k">if</span> <span class="n">cij</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span><span class="p">:</span>
                <span class="n">at</span><span class="o">.</span><span class="n">stress</span><span class="p">[:,</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">stress</span><span class="p">[:,[</span><span class="n">n</span><span class="o">.</span><span class="n">j</span> <span class="k">for</span> <span class="n">n</span> <span class="ow">in</span> <span class="n">neighb</span><span class="p">]]</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">2</span><span class="p">)</span>
                <span class="n">compute_stress_eig</span><span class="p">(</span><span class="n">i</span><span class="p">)</span>

            <span class="n">at</span><span class="o">.</span><span class="n">strain_energy_density</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">at</span><span class="o">.</span><span class="n">strain</span><span class="p">[:,</span><span class="n">i</span><span class="p">],</span> <span class="n">at</span><span class="o">.</span><span class="n">stress</span><span class="p">[:,</span><span class="n">i</span><span class="p">])</span>

    <span class="k">if</span> <span class="p">(</span><span class="n">system</span> <span class="o">==</span> <span class="s1">&#39;anatase&#39;</span><span class="p">):</span>
        <span class="nb">print</span> <span class="s1">&#39;Ti atoms computed&#39;</span><span class="p">,</span> <span class="n">n_at</span></div>


<div class="viewcode-block" id="AtomResolvedStressField"><a class="viewcode-back" href="../../elasticity.html#quippy.elasticity.AtomResolvedStressField">[docs]</a><span class="k">class</span> <span class="nc">AtomResolvedStressField</span><span class="p">(</span><span class="nb">object</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Calculator interface to :func:`elastic_fields_fortran` and :func:`elastic_fields`</span>

<span class="sd">    Computes local stresses from atom resolved strain tensor and linear elasticity.</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    bulk: Atoms object, optional</span>
<span class="sd">       If present, set `a` and ``cij`` from ``bulk.cell[0,0]`` and</span>
<span class="sd">       ``bulk.get_calculator().get_elastic_constants(bulk)``. This means</span>
<span class="sd">       `bulk` should be a relaxed cubic unit cell.</span>
<span class="sd">    a : float, optional</span>
<span class="sd">       Lattice constant</span>
<span class="sd">    cij : array_like, optional</span>
<span class="sd">       6 x 6 matrix of elastic constants :math:`C_{ij}`.</span>
<span class="sd">       Can be obtained with :meth:`.Potential.get_elastic_constants&#39;.</span>
<span class="sd">    method : str</span>
<span class="sd">       Which routine to use: one of &quot;fortran&quot; or &quot;python&quot;.</span>
<span class="sd">    **extra_args : dict</span>
<span class="sd">       Extra arguments to be passed along to python :func:`elastic_fields`, e.g.</span>
<span class="sd">       for non-cubic cells.</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span>  <span class="n">bulk</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">a</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">cij</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">method</span><span class="o">=</span><span class="s1">&#39;fortran&#39;</span><span class="p">,</span> <span class="o">**</span><span class="n">extra_args</span><span class="p">):</span>
        <span class="k">if</span> <span class="n">method</span> <span class="ow">not</span> <span class="ow">in</span> <span class="p">[</span><span class="s1">&#39;fortran&#39;</span><span class="p">,</span> <span class="s1">&#39;python&#39;</span><span class="p">]:</span>
            <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s1">&#39;method should be one of &quot;fortran&quot; or &quot;python&quot;&#39;</span><span class="p">)</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">method</span> <span class="o">=</span> <span class="n">method</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">a</span> <span class="o">=</span> <span class="n">a</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">cij</span> <span class="o">=</span> <span class="n">farray</span><span class="p">(</span><span class="n">cij</span><span class="p">)</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">extra_args</span> <span class="o">=</span> <span class="n">extra_args</span>
        <span class="k">if</span> <span class="n">bulk</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span><span class="p">:</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">a</span> <span class="o">=</span> <span class="n">bulk</span><span class="o">.</span><span class="n">cell</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span>
            <span class="n">calc</span> <span class="o">=</span> <span class="n">bulk</span><span class="o">.</span><span class="n">get_calculator</span><span class="p">()</span>
            <span class="k">if</span> <span class="n">calc</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
                <span class="k">raise</span> <span class="ne">RuntimeError</span><span class="p">(</span><span class="s1">&#39;bulk Atoms passed to AtomResolvedStressField has no calculator!&#39;</span><span class="p">)</span>
            <span class="bp">self</span><span class="o">.</span><span class="n">cij</span> <span class="o">=</span> <span class="n">farray</span><span class="p">(</span><span class="n">calc</span><span class="o">.</span><span class="n">get_elastic_constants</span><span class="p">(</span><span class="n">bulk</span><span class="p">))</span>
        

<div class="viewcode-block" id="AtomResolvedStressField.get_stresses"><a class="viewcode-back" href="../../elasticity.html#quippy.elasticity.AtomResolvedStressField.get_stresses">[docs]</a>    <span class="k">def</span> <span class="nf">get_stresses</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">atoms</span><span class="p">,</span> <span class="n">cutoff</span><span class="o">=</span><span class="mf">3.0</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        Returns local stresses on `atoms` as a ``(len(atoms), 3, 3)`` array</span>
<span class="sd">        &quot;&quot;&quot;</span>
        <span class="k">if</span> <span class="ow">not</span> <span class="nb">isinstance</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">Atoms</span><span class="p">):</span>
            <span class="n">atoms</span> <span class="o">=</span> <span class="n">Atoms</span><span class="p">(</span><span class="n">atoms</span><span class="p">)</span>
        
        <span class="n">sigma</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="nb">len</span><span class="p">(</span><span class="n">atoms</span><span class="p">),</span> <span class="mi">3</span><span class="p">,</span> <span class="mi">3</span><span class="p">))</span>
        <span class="k">if</span> <span class="bp">self</span><span class="o">.</span><span class="n">method</span> <span class="o">==</span> <span class="s1">&#39;fortran&#39;</span><span class="p">:</span>
            <span class="n">atoms</span><span class="o">.</span><span class="n">set_cutoff</span><span class="p">(</span><span class="n">cutoff</span><span class="p">)</span>
            <span class="n">atoms</span><span class="o">.</span><span class="n">calc_connect</span><span class="p">()</span>
            <span class="n">elastic_fields_fortran</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">a</span><span class="o">=</span><span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">,</span> <span class="n">cij</span><span class="o">=</span><span class="bp">self</span><span class="o">.</span><span class="n">cij</span><span class="p">)</span>
            
            <span class="n">sigma</span><span class="p">[:,</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">],</span> <span class="n">sigma</span><span class="p">[:,</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">],</span> <span class="n">sigma</span><span class="p">[:,</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">],</span> <span class="n">sigma</span><span class="p">[:,</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">],</span> <span class="n">sigma</span><span class="p">[:,</span><span class="mi">0</span><span class="p">,</span><span class="mi">2</span><span class="p">],</span> <span class="n">sigma</span><span class="p">[:,</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> \
               <span class="n">atoms</span><span class="o">.</span><span class="n">sig_xx</span><span class="p">,</span> <span class="n">atoms</span><span class="o">.</span><span class="n">sig_yy</span><span class="p">,</span> <span class="n">atoms</span><span class="o">.</span><span class="n">sig_zz</span><span class="p">,</span> <span class="n">atoms</span><span class="o">.</span><span class="n">sig_yz</span><span class="p">,</span> <span class="n">atoms</span><span class="o">.</span><span class="n">sig_xz</span><span class="p">,</span> <span class="n">atoms</span><span class="o">.</span><span class="n">sig_xy</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="n">elastic_fields</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">a</span><span class="o">=</span><span class="bp">self</span><span class="o">.</span><span class="n">a</span><span class="p">,</span> <span class="n">cij</span><span class="o">=</span><span class="bp">self</span><span class="o">.</span><span class="n">cij</span><span class="p">,</span> <span class="o">**</span><span class="bp">self</span><span class="o">.</span><span class="n">extra_args</span><span class="p">)</span>

            <span class="n">sigma</span><span class="p">[:,</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">],</span> <span class="n">sigma</span><span class="p">[:,</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">],</span> <span class="n">sigma</span><span class="p">[:,</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">],</span> <span class="n">sigma</span><span class="p">[:,</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">],</span> <span class="n">sigma</span><span class="p">[:,</span><span class="mi">0</span><span class="p">,</span><span class="mi">2</span><span class="p">],</span> <span class="n">sigma</span><span class="p">[</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">stress</span>

        <span class="c1"># Fill in symmetric components</span>
        <span class="n">sigma</span><span class="p">[:,</span><span class="mi">1</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">sigma</span><span class="p">[:,</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span>
        <span class="n">sigma</span><span class="p">[:,</span><span class="mi">2</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">sigma</span><span class="p">[:,</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span>
        <span class="n">sigma</span><span class="p">[:,</span><span class="mi">2</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">sigma</span><span class="p">[:,</span><span class="mi">0</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span>

        <span class="k">return</span> <span class="n">sigma</span></div>

<div class="viewcode-block" id="AtomResolvedStressField.get_stress"><a class="viewcode-back" href="../../elasticity.html#quippy.elasticity.AtomResolvedStressField.get_stress">[docs]</a>    <span class="k">def</span> <span class="nf">get_stress</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">atoms</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        Returns total stress on `atoms`, as a 6-element array</span>
<span class="sd">        &quot;&quot;&quot;</span>
        <span class="n">stresses</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">get_stresses</span><span class="p">(</span><span class="n">atoms</span><span class="p">)</span>
        <span class="k">return</span> <span class="n">stresses</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span></div></div>
        

<span class="k">def</span> <span class="nf">elasticity_matrix_to_tensor</span><span class="p">(</span><span class="n">C</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;Given a 6x6 elastic matrix in compressed Voigt notation, return</span>
<span class="sd">    the full 3x3x3x3 elastic constant tensor C_ijkl&quot;&quot;&quot;</span>

    <span class="n">c</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">((</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">),</span><span class="s1">&#39;d&#39;</span><span class="p">)</span>
    <span class="k">for</span> <span class="n">p</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">6</span><span class="p">):</span>
        <span class="k">for</span> <span class="n">q</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">6</span><span class="p">):</span>
            <span class="n">i</span><span class="p">,</span> <span class="n">j</span> <span class="o">=</span> <span class="n">voigt_map</span><span class="p">[</span><span class="n">p</span><span class="p">]</span>
            <span class="n">k</span><span class="p">,</span> <span class="n">l</span> <span class="o">=</span> <span class="n">voigt_map</span><span class="p">[</span><span class="n">q</span><span class="p">]</span>
            <span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">]</span> <span class="o">=</span> <span class="n">c</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">i</span><span class="p">,</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">]</span> <span class="o">=</span> <span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">l</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">=</span> <span class="n">c</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">i</span><span class="p">,</span><span class="n">l</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">=</span> <span class="n">c</span><span class="p">[</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">C</span><span class="p">[</span><span class="n">p</span><span class="p">,</span><span class="n">q</span><span class="p">]</span>

    <span class="k">return</span> <span class="n">c</span>

<span class="k">def</span> <span class="nf">elasticity_tensor_to_matrix</span><span class="p">(</span><span class="n">c</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;Given full tensor c_ijkl, return compressed Voigt form C_ij.</span>
<span class="sd">    First checks that c_ijkl obeys required symmetries.&quot;&quot;&quot;</span>

    <span class="n">tol</span> <span class="o">=</span> <span class="mf">1e-10</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
            <span class="k">for</span> <span class="n">k</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
                <span class="k">for</span> <span class="n">l</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
                    <span class="k">assert</span><span class="p">(</span><span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">]</span> <span class="o">-</span> <span class="n">c</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">i</span><span class="p">,</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">]</span> <span class="o">&lt;</span> <span class="n">tol</span><span class="p">)</span>
                    <span class="k">assert</span><span class="p">(</span><span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">]</span> <span class="o">-</span> <span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">l</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">&lt;</span> <span class="n">tol</span><span class="p">)</span>
                    <span class="k">assert</span><span class="p">(</span><span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">]</span> <span class="o">-</span> <span class="n">c</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">i</span><span class="p">,</span><span class="n">l</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">&lt;</span> <span class="n">tol</span><span class="p">)</span>
                    <span class="k">assert</span><span class="p">(</span><span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">]</span> <span class="o">-</span> <span class="n">c</span><span class="p">[</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">&lt;</span> <span class="n">tol</span><span class="p">)</span>

    <span class="n">C</span> <span class="o">=</span> <span class="n">zeros</span><span class="p">((</span><span class="mi">6</span><span class="p">,</span><span class="mi">6</span><span class="p">),</span><span class="s1">&#39;d&#39;</span><span class="p">)</span>
    <span class="k">for</span> <span class="n">p</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">6</span><span class="p">):</span>
        <span class="k">for</span> <span class="n">q</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">6</span><span class="p">):</span>
            <span class="n">i</span><span class="p">,</span> <span class="n">j</span> <span class="o">=</span> <span class="n">voigt_map</span><span class="p">[</span><span class="n">p</span><span class="p">]</span>
            <span class="n">k</span><span class="p">,</span> <span class="n">l</span> <span class="o">=</span> <span class="n">voigt_map</span><span class="p">[</span><span class="n">q</span><span class="p">]</span>
            <span class="n">C</span><span class="p">[</span><span class="n">p</span><span class="p">,</span><span class="n">q</span><span class="p">]</span> <span class="o">=</span> <span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">]</span>

    <span class="k">return</span> <span class="n">C</span>

<span class="c1"># Voigt notation: 1 = 11 (xx), 2 = 22 (yy), 3 = 33 (zz), 4 = 23 (yz), 5 = 31 (zx), 6 = 12 (xy</span>
<span class="n">voigt_map</span> <span class="o">=</span> <span class="nb">dict</span><span class="p">(</span> <span class="p">((</span><span class="mi">1</span><span class="p">,(</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">)),</span>  <span class="c1"># xx</span>
                   <span class="p">(</span><span class="mi">2</span><span class="p">,(</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">)),</span>  <span class="c1"># yy</span>
                   <span class="p">(</span><span class="mi">3</span><span class="p">,(</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">)),</span>  <span class="c1"># zz</span>
                   <span class="p">(</span><span class="mi">4</span><span class="p">,(</span><span class="mi">2</span><span class="p">,</span><span class="mi">3</span><span class="p">)),</span>  <span class="c1"># yz</span>
                   <span class="p">(</span><span class="mi">5</span><span class="p">,(</span><span class="mi">3</span><span class="p">,</span><span class="mi">1</span><span class="p">)),</span>  <span class="c1"># zx</span>
                   <span class="p">(</span><span class="mi">6</span><span class="p">,(</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">))))</span> <span class="c1"># xy</span>

<span class="k">def</span> <span class="nf">elasticity_matrix_to_tensor</span><span class="p">(</span><span class="n">C</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;Given a 6x6 elastic matrix in compressed Voigt notation, return</span>
<span class="sd">    the full 3x3x3x3 elastic constant tensor C_ijkl&quot;&quot;&quot;</span>

    <span class="n">c</span> <span class="o">=</span> <span class="n">fzeros</span><span class="p">((</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">),</span><span class="s1">&#39;d&#39;</span><span class="p">)</span>
    <span class="k">for</span> <span class="n">p</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">6</span><span class="p">):</span>
        <span class="k">for</span> <span class="n">q</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">6</span><span class="p">):</span>
            <span class="n">i</span><span class="p">,</span> <span class="n">j</span> <span class="o">=</span> <span class="n">voigt_map</span><span class="p">[</span><span class="n">p</span><span class="p">]</span>
            <span class="n">k</span><span class="p">,</span> <span class="n">l</span> <span class="o">=</span> <span class="n">voigt_map</span><span class="p">[</span><span class="n">q</span><span class="p">]</span>
            <span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">]</span> <span class="o">=</span> <span class="n">c</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">i</span><span class="p">,</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">]</span> <span class="o">=</span> <span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">l</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">=</span> <span class="n">c</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">i</span><span class="p">,</span><span class="n">l</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">=</span> <span class="n">c</span><span class="p">[</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">C</span><span class="p">[</span><span class="n">p</span><span class="p">,</span><span class="n">q</span><span class="p">]</span>

    <span class="k">return</span> <span class="n">c</span>

<span class="k">def</span> <span class="nf">elasticity_tensor_to_matrix</span><span class="p">(</span><span class="n">c</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;Given full tensor c_ijkl, return compressed Voigt form C_ij.</span>
<span class="sd">    </span>
<span class="sd">    First checks that c_ijkl obeys required symmetries.&quot;&quot;&quot;</span>

    <span class="n">tol</span> <span class="o">=</span> <span class="mf">1e-10</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
            <span class="k">for</span> <span class="n">k</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
                <span class="k">for</span> <span class="n">l</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
                    <span class="k">assert</span><span class="p">(</span><span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">]</span> <span class="o">-</span> <span class="n">c</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">i</span><span class="p">,</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">]</span> <span class="o">&lt;</span> <span class="n">tol</span><span class="p">)</span>
                    <span class="k">assert</span><span class="p">(</span><span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">]</span> <span class="o">-</span> <span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">l</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">&lt;</span> <span class="n">tol</span><span class="p">)</span>
                    <span class="k">assert</span><span class="p">(</span><span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">]</span> <span class="o">-</span> <span class="n">c</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">i</span><span class="p">,</span><span class="n">l</span><span class="p">,</span><span class="n">k</span><span class="p">]</span> <span class="o">&lt;</span> <span class="n">tol</span><span class="p">)</span>
                    <span class="k">assert</span><span class="p">(</span><span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">]</span> <span class="o">-</span> <span class="n">c</span><span class="p">[</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">,</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">&lt;</span> <span class="n">tol</span><span class="p">)</span>

    <span class="n">C</span> <span class="o">=</span> <span class="n">fzeros</span><span class="p">((</span><span class="mi">6</span><span class="p">,</span><span class="mi">6</span><span class="p">),</span><span class="s1">&#39;d&#39;</span><span class="p">)</span>
    <span class="k">for</span> <span class="n">p</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">6</span><span class="p">):</span>
        <span class="k">for</span> <span class="n">q</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">6</span><span class="p">):</span>
            <span class="n">i</span><span class="p">,</span> <span class="n">j</span> <span class="o">=</span> <span class="n">voigt_map</span><span class="p">[</span><span class="n">p</span><span class="p">]</span>
            <span class="n">k</span><span class="p">,</span> <span class="n">l</span> <span class="o">=</span> <span class="n">voigt_map</span><span class="p">[</span><span class="n">q</span><span class="p">]</span>
            <span class="n">C</span><span class="p">[</span><span class="n">p</span><span class="p">,</span><span class="n">q</span><span class="p">]</span> <span class="o">=</span> <span class="n">c</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">]</span>

    <span class="k">return</span> <span class="n">C</span>

<div class="viewcode-block" id="transform_elasticity"><a class="viewcode-back" href="../../elasticity.html#quippy.elasticity.transform_elasticity">[docs]</a><span class="k">def</span> <span class="nf">transform_elasticity</span><span class="p">(</span><span class="n">c</span><span class="p">,</span> <span class="n">R</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;Transform c as a rank-4 tensor by the rotation matrix R.</span>

<span class="sd">    Returns the new representation c&#39;. If c is a 6x6 matrix it is first</span>
<span class="sd">    converted to 3x3x3x3 form, and then converted back after the</span>
<span class="sd">    transformation.&quot;&quot;&quot;</span>

    <span class="n">made_tensor</span> <span class="o">=</span> <span class="kc">False</span>
    <span class="k">if</span> <span class="n">c</span><span class="o">.</span><span class="n">shape</span> <span class="o">==</span> <span class="p">(</span><span class="mi">6</span><span class="p">,</span><span class="mi">6</span><span class="p">):</span>
        <span class="n">made_tensor</span> <span class="o">=</span> <span class="kc">True</span>
        <span class="n">c</span> <span class="o">=</span> <span class="n">elasticity_matrix_to_tensor</span><span class="p">(</span><span class="n">c</span><span class="p">)</span>
    <span class="k">elif</span> <span class="n">c</span><span class="o">.</span><span class="n">shape</span> <span class="o">==</span> <span class="p">(</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">):</span>
        <span class="k">pass</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s1">&#39;Tensor should either be 3x3x3x3 or 6x6 matrix&#39;</span><span class="p">)</span>

    <span class="n">cp</span> <span class="o">=</span> <span class="n">fzeros</span><span class="p">((</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">),</span><span class="s1">&#39;d&#39;</span><span class="p">)</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
            <span class="k">for</span> <span class="n">k</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
                <span class="k">for</span> <span class="n">l</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
                    <span class="k">for</span> <span class="n">p</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
                        <span class="k">for</span> <span class="n">q</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
                            <span class="k">for</span> <span class="n">r</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
                                <span class="k">for</span> <span class="n">s</span> <span class="ow">in</span> <span class="n">frange</span><span class="p">(</span><span class="mi">3</span><span class="p">):</span>
                                    <span class="n">cp</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">,</span><span class="n">k</span><span class="p">,</span><span class="n">l</span><span class="p">]</span> <span class="o">+=</span> <span class="n">R</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">p</span><span class="p">]</span><span class="o">*</span><span class="n">R</span><span class="p">[</span><span class="n">j</span><span class="p">,</span><span class="n">q</span><span class="p">]</span><span class="o">*</span><span class="n">R</span><span class="p">[</span><span class="n">k</span><span class="p">,</span><span class="n">r</span><span class="p">]</span><span class="o">*</span><span class="n">R</span><span class="p">[</span><span class="n">l</span><span class="p">,</span><span class="n">s</span><span class="p">]</span><span class="o">*</span><span class="n">c</span><span class="p">[</span><span class="n">p</span><span class="p">,</span><span class="n">q</span><span class="p">,</span><span class="n">r</span><span class="p">,</span><span class="n">s</span><span class="p">]</span>

    <span class="k">if</span> <span class="n">made_tensor</span><span class="p">:</span>
        <span class="k">return</span> <span class="n">elasticity_tensor_to_matrix</span><span class="p">(</span><span class="n">cp</span><span class="p">)</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="k">return</span> <span class="n">cp</span></div>


<div class="viewcode-block" id="rayleigh_wave_speed"><a class="viewcode-back" href="../../elasticity.html#quippy.elasticity.rayleigh_wave_speed">[docs]</a><span class="k">def</span> <span class="nf">rayleigh_wave_speed</span><span class="p">(</span><span class="n">C</span><span class="p">,</span> <span class="n">rho</span><span class="p">,</span> <span class="n">a</span><span class="o">=</span><span class="mf">4000.</span><span class="p">,</span> <span class="n">b</span><span class="o">=</span><span class="mf">6000.</span><span class="p">,</span> <span class="n">isotropic</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Rayleigh wave speed in a crystal.</span>

<span class="sd">    Returns triplet ``(vs, vp, c_R)`` in m/s, where `vs` is the transverse</span>
<span class="sd">    wave speed, `vp` the longitudinal wave speed and `c_R` the Rayleigh</span>
<span class="sd">    shear wave speed. </span>

<span class="sd">    For the anisotropic case (default), formula is Darinskii,</span>
<span class="sd">    A. (1997).  On the theory of leaky waves in crystals.  `Wave</span>
<span class="sd">    Motion, 25(1),</span>
<span class="sd">    35-49. &lt;http://dx.doi.org/10.1016/S0165-2125(96)00031-5&gt;`_. If</span>
<span class="sd">    `isostropic` is True, formula is from `this page</span>
<span class="sd">    &lt;http://sepwww.stanford.edu/data/media/public/docs/sep124/jim1/paper_html/node5.html&gt;`_</span>

<span class="sd">    `C` is the 6x6 matrix of elastic contstant, rotated to reference</span>
<span class="sd">    the frame of sample, and should be given in units of GPa.  The</span>
<span class="sd">    Rayleight speed returned is along the first (x) axis.</span>

<span class="sd">    `rho` is the density in g/cm^3.</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="k">if</span> <span class="s1">&#39;scipy&#39;</span> <span class="ow">not</span> <span class="ow">in</span> <span class="n">available_modules</span><span class="p">:</span>
        <span class="k">raise</span> <span class="ne">RuntimeError</span><span class="p">(</span><span class="s1">&#39;scipy is needed for rayleigh_wave_speed()&#39;</span><span class="p">)</span>

    <span class="kn">from</span> <span class="nn">scipy.optimize</span> <span class="k">import</span> <span class="n">bisect</span>

    <span class="n">C</span> <span class="o">=</span> <span class="n">C</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>
    <span class="n">C</span> <span class="o">*=</span> <span class="mf">1e9</span> <span class="c1"># convert GPa -&gt; Pa</span>
    <span class="n">rho</span> <span class="o">*=</span> <span class="mf">1e3</span> <span class="c1"># convert g/cm^3 -&gt; kg/m^3</span>

    <span class="k">if</span> <span class="n">isotropic</span><span class="p">:</span>
        <span class="n">vp</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">C</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span><span class="o">/</span><span class="n">rho</span><span class="p">)</span>
        <span class="n">vs</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">C</span><span class="p">[</span><span class="mi">4</span><span class="p">,</span><span class="mi">4</span><span class="p">]</span><span class="o">/</span><span class="n">rho</span><span class="p">)</span>
        <span class="n">f</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">v</span><span class="p">:</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="p">(</span><span class="n">v</span><span class="o">/</span><span class="n">vs</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="p">(</span><span class="n">v</span><span class="o">/</span><span class="n">vp</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span> <span class="o">-</span>
                       <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">v</span><span class="o">**</span><span class="mi">2</span><span class="o">/</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">vs</span><span class="o">**</span><span class="mi">2</span><span class="p">)))</span><span class="o">**</span><span class="mi">4</span><span class="p">)</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">vp</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">C</span><span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">]</span><span class="o">/</span><span class="n">rho</span><span class="p">)</span>
        <span class="n">vs</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">C</span><span class="p">[</span><span class="mi">4</span><span class="p">,</span><span class="mi">4</span><span class="p">]</span><span class="o">/</span><span class="n">rho</span><span class="p">)</span>
        <span class="n">f</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="p">((</span><span class="n">x</span><span class="o">/</span><span class="n">vp</span><span class="p">)</span><span class="o">**</span><span class="mi">4</span><span class="o">*</span><span class="p">(</span> <span class="p">(</span><span class="n">C</span><span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">]</span><span class="o">/</span><span class="n">C</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">])</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="p">(</span><span class="n">x</span><span class="o">/</span><span class="n">vp</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">))</span> <span class="o">-</span>
                       <span class="p">((</span><span class="mi">1</span><span class="o">-</span><span class="p">(</span> <span class="p">(</span><span class="n">C</span><span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">C</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span><span class="o">*</span><span class="n">C</span><span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">])</span> <span class="p">)</span> <span class="o">-</span>
                       <span class="p">(</span><span class="n">x</span><span class="o">/</span><span class="n">vp</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">*</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="p">(</span><span class="n">x</span><span class="o">/</span><span class="n">vs</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">))</span>
    <span class="n">c_R</span> <span class="o">=</span> <span class="n">bisect</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="n">vs</span><span class="o">/</span><span class="mf">2.</span><span class="p">,</span> <span class="n">vs</span><span class="p">)</span>
    <span class="k">return</span> <span class="n">vs</span><span class="p">,</span> <span class="n">vp</span><span class="p">,</span> <span class="n">c_R</span></div>

</pre></div>

           </div>
           
          </div>
          <footer>
  

  <hr/>

  <div role="contentinfo">
    <p>
        &copy; Copyright 2008-2016, James Kermode

    </p>
  </div>
  Built with <a href="http://sphinx-doc.org/">Sphinx</a> using a <a href="https://github.com/rtfd/sphinx_rtd_theme">theme</a> provided by <a href="https://readthedocs.org">Read the Docs</a>. 

</footer>

        </div>
      </div>

    </section>

  </div>
  


  <script type="text/javascript">
      jQuery(function () {
          SphinxRtdTheme.Navigation.enable(true);
      });
  </script>

  
  
    
   

</body>
</html>